####################################
# Main effect RDD
####################################

rm(list=ls())

library(gridExtra)
library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)

################
# Prepare data 
################

# read data
d = read.dta13("~/Dropbox/Anti-american attitudes/01_data/clean_data/trump_election_data_clean_final.dta")   
names(d)

# baseline support
mean(d$us_trust_bin[d$treatment==0])
round(tapply(d$us_trust_bin[d$treatment==0],d$country[d$treatment==0],mean),2)

######################
# US trust binary 
######################

# bandwidths
h_us_trust_bin = rdbwselect(d$us_trust_bin,d$score)$bws[1]
h_us_trust_bin

# optimal 
kweights = kernelwts(d$score, 0, bw = h_us_trust_bin, kernel = "triangular")
coef(summary(lm(d$us_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)

tableresults = cbind(results$days,results$n,results$pe,results$se,results$pv)
tableresults

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

######################
# Missingness
######################

# bandwidths
h_miss = rdbwselect(d$missigness,d$score)$bws[1]
h_miss

# optimal 
kweights = kernelwts(d$score, 0, bw = h_miss, kernel = "triangular")
coef(summary(lm(d$missigness ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$missigness ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$missigness ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$missigness ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$missigness))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_miss))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)

tableresults = cbind(results$days,results$n,results$pe,results$se,results$pv)
tableresults

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

######################
# US trust continuous
######################

# subset
table(d$us_trust_con, exclude = NULL)
d2 = d[!is.na(d$us_trust_con),]

# bandwidth
h_us_trust_con = rdbwselect(d2$us_trust_con,d2$score)$bws[1]
h_us_trust_con

# optimal
kweights = kernelwts(d2$score, 0, bw = h_us_trust_con, kernel = "triangular")
coef(summary(lm(d2$us_trust_con ~ d2$treatment + d2$score + d2$treatment*d2$score + as.factor(d2$pais), weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d2$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d2$us_trust_con ~ d2$treatment + d2$score + d2$treatment*d2$score + as.factor(d2$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d2$us_trust_con ~ d2$treatment + d2$score + d2$treatment*d2$score + as.factor(d2$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d2$us_trust_con ~ d2$treatment + d2$score + d2$treatment*d2$score + as.factor(d2$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d2$us_trust_con))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_con))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
round(results$pv,3)

tableresults = cbind(results$days,results$n,results$pe,results$se,results$pv)
tableresults

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4)) + coord_cartesian(ylim=c(-0.4,0.4)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

######################
# AF trust binary 
######################

# bandwidth
h_af_trust_bin = rdbwselect(d$af_trust_bin,d$score)$bws[1]
h_af_trust_bin

# optimal 
kweights = kernelwts(d$score, 0, bw = h_af_trust_bin, kernel = "triangular")
coef(summary(lm(d$af_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$af_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$af_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$af_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$af_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_af_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
mean(results$pv,3)

tableresults = cbind(results$days,results$n,results$pe,results$se,results$pv)
tableresults

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 

#############################
# Municipality trust binary 
#############################

# bandwidth
h_muni_trust_bin = rdbwselect(d$muni_trust_bin,d$score)$bws[1]
h_muni_trust_bin

# optimal 
kweights = kernelwts(d$score, 0, bw = h_muni_trust_bin, kernel = "triangular")
coef(summary(lm(d$muni_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2,]

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$muni_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$muni_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$muni_trust_bin ~ d$treatment + d$score + d$treatment*d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$muni_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_muni_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)
mean(results$pv,3)

tableresults = cbind(results$days,results$n,results$pe,results$se,results$pv)
tableresults

ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") 
